Multiscale Analysis of the Stress State in a Granular Slope in Transition to Failure 
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By means of contact dynamics simulations, we analyze the stress state in a granular bed slowly 
tilted towards its angle of repose. An increasingly large number of grains are overloaded in the sense 
that they are found to carry a stress ratio above the Coulomb yield threshold of the whole packing. 
Using this property, we introduce a coarse-graining length scale at which all stress ratios are below 
the packing yield threshold. We show that this length increases with the slope angle and jumps 
to a length comparable to the depth of the granular bed at an angle below the angle of repose. 
This transition coincides with the onset of dilatation in the packing. We map this transition into 
a percolation transition of the overloaded grains, and we argue that in the presence of long-range 
correlations above the transition angle, the granular slope is metastable. 
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I. INTRODUCTION 

The science of granular materials was initiated by 
Coulomb's analysis of the equilibrium and failure of a 
granular talus The well-known Coulomb's failure 

criterion was later incorporated in the framework of 
a rigid-plastic behaviour based on experimental testing 
of granular samples with homogeneous boundary condi- 
tions [3-0-tH- Two centuries after Coulomb, the slope 
failure phenomena continue to interest scientists from 
various fields with evident applications to geologicalpro- 
cesses and industrial handling of granular materials 0.0- 
The main reason is that the phenomena involved in the 
evolution of a granular slope are richer than what might 
be expected from a mean macroscopic analysis 0, H H ■ 
On the other hand, the mechanisms that lead to slope 
failure have not yet been well understood from a grain 
scale standpoint [H El El El- 

New investigation tools, such as fine imaging tech- 
niques and discrete numerical simulations, have shown 
that granular media are very inhomogeneous at the grain 
scale, and the microstructure, i.e. the organization of the 
grains and their contacts in space, can evolve in many dif- 
ferent ways in response to external loading and boundary 
conditions El El El El E3 Often, large fluctua- 
tions are observed in the course of shearing and from one 
situation to another J2jj. Several observations suggest 
that surface failure may occur at slope angles well below 
the angle of repose [21|, and metastable states exist in 
the vicinity of the angle of repose . 

In this context, a closer look at the microstructure and 
spatio-temporal scales governing the behaviour of a gran- 
ular slope can not be avoided. The query is which inter- 
nal variables or order parameters represent the evolution 
of a granular slope towards surface failure [2^ • Even un- 
der "quasi-static" conditions, the grains in a cohesionless 
granular medium exhibit a high degree of mobility. Both 
dynamical instabilities and collective rearrangement phe- 
nomena occur frequently in response to slightest load in- 
crements [la, [23, |23 . This observation shows that, even 



in a granular medium far from macroscopic failure, the 
failure conditions are often fulfilled locally. Such effects 
may be observed by looking at grain displacements or 
contact forces at different scales. 

In this paper, we focus on the scaling of local stresses 
in a two-dimensional granular bed simulated by the Con- 
tact Dynamics method. The bed is tilted slowly towards 
its angle of repose C at which slope failure occurs. This 
slow gravity loading gives rise to grain rearrangements 
that harden the bed. Without such a hardening (or plas- 
tification) process, the bed can not reach its angle of 
repose. For example, a bed prepared initially by pour- 
ring the grains onto a rough horizontal plane, when tilted 
suddenly to a finite slope 8, will fail immediately. Hence, 
the evolution of local stress states as a function of load- 
ing is important for understanding how the failure limit 
is reached and in which respects it is controlled by the 
details of the microstructure. 

For the analysis of stresses at different scales, we 
need a quantity that plays the same role as the Cauchy 
stress tensor in continuous media. It can be shown that 
the concept of "internal moment tensor" , introduced by 
Moreau [23j, generalizes consistently the Cauchy stress 
tensor to discrete media and, what is more, its mechani- 
cal content with respect to Newton's equations of motion 
remains the same whether applied to a single grain or to 
a collection of grains inside a control volume. 

In the following, we first introduce the numerical pro- 
cedures and the concept of stress tensor in terms of inter- 
nal moments which will be used throughout the paper. 
Then, we analyze the data from a global standpoint, i.e. 
by considering the evolution of average stresses in the 
granular bed as a whole. We discuss also a theoreti- 
cal evaluation of the stresses in the presence of walls in 
comparison to a bed with periodic boundary conditions. 
Then, we analyze local and coarse-grained stresses as a 
function of the tilt angle. We will conclude with a sum- 
mary of our main results and a discussion about their 
physical interpretation. 
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II. NUMERICAL PROCEDURES 

A. Simulation method 

For a discrete simulation of the dynamics of a collec- 
tion of cohesionless rigid grains, two strategies can be 
adopted, differring mainly in the implementation of con- 
tact repulsion and friction. In the popular Molecular 
Dynamics (MD) method a repulsive potential is intro- 
duced as a function of contact interpenetration and the 
friction force obeys an elastic or viscous law up to a 
Coulomb threshold pBl l26l . The equations of mo- 
tion are explicitely integrated according to different clas- 
sical schemes. As an interesting alternative to the MD 
method, the Contact Dynamics (CD) approach deals di- 
rectly with infinitely stiff contact laws or, more generally 
"nonsmooth" laws [H |H E<J|3(J El H3 Hence, in con- 
trast to the MD method, no elastic stiffness or viscous 
regularization of the Coulomb friction law are to be in- 
troduced at the contact level. This allows for larger time 
steps and thus for a substential reduction of computa- 
tional time. At the same time, no purely computational 
parameters are introduced, the only parameters being 
the intergrain coefficient of friction \x and the normal and 
tangential coefficients of restitution that account for dis- 
sipation in the advent of collisions between the grains. 
In application to quasi-static deformations of a granu- 
lar packing, one expects similar behaviours from both 
methods as far as large time scales (beyond the elastic 
response time) and plastic deformations are concerned. 
In the present work, we applied the CD method. 

B. System characteristics 

We consider circular grains in two dimensions, 
with diameters uniformly distributed in the interval 
[D m in, D max ] with D max /D min = 1.5. This slight poly- 
dispersity reduces long-range cristal-like ordering in the 
packing. The grains interact only through frictional co- 
hesionless contacts with a coefficient of friction /i = 0.5. 
The same value is used for grain-wall contacts. More- 
over, we assume zero coefficients of restitution between 
grains. Slightly larger coefficients of restitution do not 
influence the mechanical behaviour in a dense granular 
packing where multiple contacts dissipate efficiently the 
kinetic energy. 

The granular beds were prepared by random deposition 
of grains on a horizontal plane in the gravity field. The 
plane was made rough by sticking grains of diameter D, 
D being the mean diameter of the grains in the packing. 
Two different boundary conditions were implemented: 
wall boundary conditions (WBC) where the bed is con- 
fined between two vertical walls, and periodic boundary 
conditions (PBC) in the horizontal direction. We pre- 
pared different beds having all the same depth H ~ 40-D, 
but different lengths L. In the following, we will ana- 
lyze WBC systems with L = 100£>, 150£>, 200D, 250D 




FIG. 1: Schema of the simulation: a granular bed with wall 
boundary conditions is tilted in the gravity field to bring the 
slope of the free surface 9 from originally up to the angle of 
avalanche 6 C . 



and 300£>, corresponding to 4000, 6000, 8000, 10000 and 
12000 grains respectively, as well as a PBC system com- 
posed of 8000 grains with L — 200, where L refers in that 
case to the length of the simulation cell. 

The general features of the packings are the same inde- 
pendently of the boundary conditions. They have a solid 
fraction p ~ 0.8 and a coordination number z ~ 3.5, 
corresponding to a random close packing. Despite the 
polydispersity introduced in the grain size distribution, 
we observe privileged directions of contact normals at 0° , 
60° and 120° with respect to the horizontal direction. A 
signature of this ordering appears in the distribution of 
local stresses, as we shall see below. 

The granular beds are tilted at a constant rotation rate 
u) = l°s _1 in the gravity field. The slope of the free sur- 
face 9 increases monotonically from 9 = to the angle 
of repose 9 = 9 C (Figure^. At this point, the stability 
limit of the packing is reached and a surface flow is trig- 
gered. Analyzing a set of 50 independent simulations, we 
checked that both the average and local behaviours that 
will be discussed in this paper are highly reproducible. 



III. THE STRESS TENSOR 

The dynamics of a granular system is naturally de- 
scribed in terms of grain degrees of freedom (velocities) 
and contact actions, including normal and tangential 
forces as well as contact torques (in a cohesive granu- 
lar medium). This vectorial description is, however, un- 
suitable in a macroscopic formulation of the rheological 
behavior where the material has to be described as an 
effective continuous medium where the stress and strain 
variables are stress tensors <r and strain tensors e de- 
fined over representative volumes (in contrast to forces 
and velocities which act over points). 

In principle, it is not difficult to evaluate the stress 
components in a control volume by simply calculating 
the surface density of forces applied by the grains situ- 
ated at one side of a surface (a line in 2D) to the grains 
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FIG. 2: The stress tensor in the volume bounded by the 
dashed line can be evaluated as the surface density of the 
forces exerted by the external grains on the grains located 
inside the volume. 



situated on the other side (Figure^. In this sense, the 
Cauchy stress tensor is as well defined as in other ma- 
terials (with or without a granular structure). In the 
case of dilute granular materials, the convection of grain 
momenta is the main mechanism of stress transmission 
as in a classical gaz. Here, we are concerned only with 
quasi-static deformations of a granular material for which 
contact actions play the main role in stress transmission. 
Using the above operational definition, it is possible to 
extract an expression for the stress components <xy in 
terms of contact forces {/"} and the vectors £ a joining 
the particle centers: 



(i) 



where V is the control volume and a denotes the contacts 
in V. 

Several authors have proposed different methods to 
demonstrate this or similar relations 

SHE!- such 

expressions are well-defined if evaluated over a large vol- 
ume containing many contacts. 

There is, of course, no reason for not using the ex- 
pression ^ at smaller scales down to a single grain if the 
volume V is properly adjusted. However, strictly speak- 
ing, the result will not be a stress tensor in the opera- 
tional sense defined above. This means that the problem 
should be posed in reverse order: Is there a quantity 
whose physical content remains the same whether ap- 
plied to a single grain or to a collection of grains inside a 
control volume and that tends to the Cauchy stress ten- 
sor at large scales? Moreau showed that the concept of 
internal moment tensor fulfills these conditions |23| . For 
clarity, we briefly introduce this concept below. 

In the framework of the virtual power formalism, a 
force (in the general sense) experienced by a bounded 



portion S of a material system is defined through the ex- 
pression of the power V that it develops when subjected 
to a virtual velocity field v(r). Let v(r) be an affine field, 



Vi(r) = Vi(0) + bi^rj, 



(2) 



where we assume Eintein's summation rule over sub- 
scripts. By definition, the power Vint (v) of internal forces 
is linear in v . This means that there exist R and M such 
that 



RMO) 



Mijbij 



(3) 



In the particular case of a rigid body motion, b is anti- 
symmetric (bij = — bji) and Vint = by virtue of New- 
ton's third law. This implies that R = and M is a 
symmetric tensor of rank 2 and independent of the choice 
of the reference frame. Following Moreau, we will refer 
to M as the internal moment tensor of the system [23j . 

By definition, the internal moment tensor makes sense 
at all scales. In particular, we may evaluate the inter- 
nal moment tensor of a grain within a granular system. 
The simplest example occurs when the system is in static 
equilibrium. In this case, the total power V = Vi„t+V e xt, 
where V ex t is the power associated with external forces, is 
zero independently of the choice of the virtual velocities. 
If the only forces f a acting on a grain p are those exerted 
at its contact points r a by neighbouring grains, then the 
internal power is Vi n t (p) = -V ex t(p) = ~ J2 a e P v i( ra )fi 
(Figure[3J). Identifying this with the general expression^ 
of the internal power, we get My(p) = — J2aep r ? ■fj' ■ 
If the condition of equilibrium does not apply, the total 
virtual power is given by J j(r)dm, where 7(7*) is the 
acceleration field and dm denotes the mass measure. 

In the case of circular grains with moment of inertia 
/ about the grains centers, the general expression of the 
internal moment tensor of a grain p becomes pjjj 



1 



-Ilu Si 



(4) 



where u is the rotation velocity and Sij is the Kronecker 
symbol. It can be shown that the expression 0] holds also 
in the presence of bulk forces (gravity) acting at grain 
centers if the origin of coordinates for each grain is placed 
at its center. 

The internal moment M(pi UP2) of two grains p\ and 
P2 sharing a contact is the sum of their respective inter- 
nal moments M(pi) and (p 2 ) because opposite reaction 
forces of equal magnitude act on the two grains at the 
same contact point. This additive property implies that 
the total internal moment M(S) of a system S is simply 
the sum of the internal moments of all grains included 
in S. On the other hand, if the number of grains in S is 
sufficiently large, it makes sense to evaluate the Cauchy 
stress tensor <j for S. Assuming the same test field as El 
the corresponding internal power by definition of <r is 



Vint = I OijdiVjdV. 



(5) 
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average stress P — o\ + o%: 



FIG. 3: Force f a applied on a grain p at the contact point 
r a by a neighbouring grain. 



Then, according to EI we have 



My (S) 



dijdV = {a tj )V. 



(6) 



This shows that the internal moment tensor of S per unit 
volume (M/V) tends to the average Cauchy stress tensor 
{<Jij ) at larger scales or for an increasing number of grains 
contained in S. 

We see that the internal moment tensor has all the 
required properties for a scaling analysis of stress trans- 
mission in granular media. Conceptually, the internal 
moment tensor per unit volume in a discrete system plays 
the same role as the Cauchy stress tensor in a continuous 
medium. For this reason, we will refer to the internal 
moments of the grains as grain stresses. In the following, 
we will be more specifically interested in stress ratios, i.e. 
the ratio of the deviatoric part of grain stresses normal- 
ized by the corresponding spherical part. 



IV. THE PACKING STRESS 

In this section, we study the average (or total) stress of 
the granular bed as a function of the tilt angle 9. In the 
absence of plastification, i.e. for a rigid- plastic behaviour 
without rearrangements, the average stress tensor er fol- 
lows directly the rotation of the bed with respect to the 
direction of gravity (and the volume of the bed remains 
constant) 22]. Here, we check whether the simulations 
yield a picture close to this prediction in spite of rear- 
rangement phenomena as the bed is tilted. The packing 
stress tensor is evaluated for the whole packing by sim- 
ply adding up the grain stresses and dividing by the total 
volume of the packing. We consider more specifically the 
ratio of the normal component cjn of the stress tensor 
along the direction of the free surface to the tangential 
component &t- Alternatively, we evaluate the stress ra- 
tio r defined as the stress deviator Q = o\ — 02, where 
<7i and (T2 are the principal values, normalized by the 



Q 

r = — . 

p 



(7) 



The Coulomb criterion implies T = sin# C) where 9 C is the 
angle of repose, at incipient failure. 

Figure ^displays the evolution of the normalized shear 
stress (Jt/&n as a function of tan 9 for all samples, as well 
as the analytical fit (Jt/&n = tan 9 for the PBC system 
(see below) . The shear stress ot / &n along the bed is ini- 
tially zero for all samples. It increases linearly with tan# 
but with different slopes for all samples. For all WBC 
systems, the slope is below that of the PBC system, but 
the slope increases with the length L of the granular bed. 
The angle of repose 9 C , where surface failure is initiated, 
varies between ss 19° and « 21° for the WBC and the 
PBC systems. 

The same trends are observed in Figure \E\ which shows 
the stress ratio T as a function of sin 9 for all samples. 
The stress ratio T is ~ 0.08 at 9 = 0. This corresponds to 
a slightly anisotropic stress state of the packings in the 
initial state. 

In order to evaluate the influence of the walls, let us 
consider Euler's equations for a medium in static equilib- 
rium [36( . In a reference frame attached to the box with 
its origin located at the left uppermost point of the bed 
and with its axes x and y oriented along and perpendic- 
ular to the bed (Figure QJ, respectively, we have 



d x a xx + dy<r yx = wsin9, 
d x a xy + d y a yy = wcos9, 



(8) 
(9) 



where w is the specific weight of the medium. The 
stresses are zero at the free surface. Since the stresses 
are evaluated for the whole bed, we have ctjv = (&yy) 
and ot = (fyx), where (...) denotes averaging over the 
whole bed. 

In the PBC case, the bed is translationally invariant 
along the free surface, so that d x a xx = and d x a xy = 0. 
Then, Euler's equations yield a yx — wsin^ y and a yy = 
w cos 9 y, so that 



(T'r 

<Jn 



= tan ( 



(10) 



This classical result fits excellently our data in the PBC 
case (Figure0J). In particular, at failure <7t/<tat = tan# c . 
From Eq. njit is easy to deduce also T 



T = sh^ 



(11) 



This form fits well the data in the PBC case as shown in 
Figure |3] 

In the WBC case, the translational invariance is broken 
in the presence of the walls. In order to close Euler's 
equations, we need to make an assumption about the 
stress state. Such an assumption represents either the 
outcome of a stress-stress relation, which is a standard 
method to close balance equations, or the history of the 
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FIG. 4: Evolution of the normalized shear stress (7t/&n as a 
function of tan# for all WBC and PBC systems. The dotted 
line corresponds to the analytical fit gt/&n = tan((9). 
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FIG. 5: Evolution of the stress ratio T as a function of sin# 
for all WBC and PBC systems. The dotted line corresponds 
to the analytical fit V — sin(#). 



packing, i.e. the plastification of the bed as it is rotated 
from its initial state. Our data show that the main effect 
of the walls is to introduce a globally nonzero gradient 
of the xx component along the bed (as a function of x). 
Figure shows that this effect is progressively localized 
in the vicinity of the walls as 9 is increased. At the same 
time, the other components may be considered as nearly 
independent of x up to wall effects. 

In order to capture the influence of a stress gradient 
along the bed in an analytical approach, let us simply 
assume that a xx has a constant gradient along the bed 
and d x a xy = 0. The latter together with equation |5| 
implies that a yy — w cos 9 y. The influence of higher 
order gradients may be evaluated as well. But, we are 
concerned here only with the simplest level of description. 
Now, let us use two coefficients ko and Ul to specify the 
boundary values of a xx : 



Cxx(% = 0) = kgw sin 9, 
a xx (x = L) = kLWsm.9. 



(12) 
(13) 




FIG. 6: Evolution of the a xx component of the stress tensor 
computed over successive sections of a WBC granular bed 
(L — 100D) as a function of the distance x of the section to 
the left wall, and for 3 values of the slope angle 9. 



case, we have ko — k^ — 1, and their values in the pres- 
ence of the walls should reflect thus small perturbation 
to translational invariance. 

We naturally expect that the normal stress u xx (x = L) 
on the lower wall is larger than the normal stress a xx [x = 
0) on the upper wall when 9 > 0, so that Ul > ko- Given 
these boundary conditions, and since a xx is assumed to 
have a constant gradient, we have 



&xx = {ko + j^( k L - fco)} wsin9 y. 
Then, from equation [S] we get 



<7t 
ctjv 



= tan ( 



1 - 



k L - k H 



(14) 



(15) 



These two coefficients might depend on 9. In the PBC 



This equation, with basically one fitting parameter 
— ko, captures all the features observed in Figure 0] 
In particular, since k^ — ko > 0, the ratio increases with 
L and it tends to tan 9 as H/L — » oo. Figure shows 
the variation of {ctt/vn)/ tan (9 for 9 — 9 C as a function 
of L from the numerical data, together with the ana- 
lytical form Eq. ^3 which fits correctly the data with 
fcx — fco ~ 2.5. 

In the forthcoming sections, we focus on the PBC sys- 
tem, though all the observed features are basically the 
same in WBC systems. 



V. LOCAL STRESSES 

A. Grain stress distributions 

According to the definition of grain stresses (equa- 
tion 2| , one can attribute to each grain in the packing 
a stress deviator q and a mean stress p. Figure |S] displays 
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FIG. 7: Ratio {ot/on)/ tan 8 for 8 — 8 C as a function of L for 
the WBC system (black circles). The dotted line represents 
the analytical form given by Ea. 1151 with fcx — fco — 2.5. 
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FIG. 8: Probability density function (pdf) P(pi) of the grain- 
stress ratios 7 = q/p for 8 — 0, 6 — 10° and 6 — 15° in the 
PBC system. The vertical line shows the value of the critical 
stress ratio T c of the packing. 



the probability density function (pdf) P("j) of the grain 
stress ratios 7 = q/p for = 0, 6 = 10° and 9 = 15° 
in the PBC system. These distributions are quite wide 
and cover all possible values of 7 in the range [0,1]. In 
particular, even at low slope angles, a large fraction of 
the grains has a stress ratio 7 above the critical stress 
ratio r c of the packing (represented by a vertical line in 
Figure [SJ. The initial fraction (at 9 = 0) of these "over- 
loaded" grains is about 30%. Indeed, the grain failure 
thresholds (the largest permissible grain stress ratios) are 
controlled by the immediate environment of each grain 
(orientations of contact neighbours) and hence they can 
overpass the threshold r c for the packing as a whole. This 
fraction grows with 9 and eventually reaches ~ 60% at 9 C , 
as shown in Figure El This suggests that the proportion 
of overloaded grains is a good indicator of the "distance" 
to slope failure (an order parameter, as discussed below). 

Figure ITUI displays the maps of grain stress ratios 7 at 
different stages of the evolution using a gray level scale. 
Two cutoff values are (arbitrarily) introduced in order 
to improve the contrast: The grains with 7 < 0.05 ap- 
pear in white, whereas the grains with 7 > 2T C appear in 
black. Intermediate values of 7 are encoded by a linear 
gray level scale. We can see on successive maps an in- 
creasing number of highly loaded grains that remain so in 
spite of recurrent rearrangment phenomena in the pack- 
ing. Furthermore we observe a remarkable organization 
of overloaded grains in the packing with a trend to form 
clusters. The map reveals also privileged directions (at 
60° and 120° with the horizontal direction) correspond- 
ing with higher concentrations of overloaded grains. 




e (deg) 

FIG. 9: Evolution of the fraction p of overloaded grains in 
the PBC system as a function of the slope angle 8. 

B. A percolation process 

We characterize here the clustering of overloaded 
grains by analyzing their neighbourhoods. The critical 
stress ratio T c is a large-scale property whereas the grain 
stress ratios 7 describe the grain scale stress state of the 
packing. Since we have 7 > T c for overloaded grains, 
we expect that the average stress ratio over a larger and 
larger neighbourhood of an overloaded grain will decline 
and eventually go below T c . The extent of the neighbour- 
hood for which 7 = T c is thus an intermediate length 
scale (between the grain size and the system size) that 
reflects the spatial correlations of overloaded grains. 

In order to evaluate these correlations from the data, 
we consider circular neighbourhoods of diameter £ cen- 
tered on overloaded grains (Figure fTTfl . For each grain, 
we calculate the stress ratio over its circular neigh- 
bourhood as a function of I. Thereby we determine, for 
each grain and at different slope angles 9, the length 
£ c for which 7^ = T c . Here, we study the evolution of 
the maximum value £™ ax and the mean value £™ ean = 
(1/N ) J2n ^ c as a function of 9. 

The plots of £™ ean and £™ ax as a function of 9 are 
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FIG. 10: Maps of the grains stress ratios 7 for successive 
values of 9. From top to bottom, 6 = 0, 9 — 5°, 9 — 10°, 
9 — 15° and 9 — 9 C . In black are represented the grains such 
that 7 > 2r c , while the white grains corresponds to 7 < 0.05; 
intermediate values are encoded by a linear gray scale. 

displayed in Figure El The mean radius £™ ean increases 
very slowly from AD to 5D. Then, from 9 = 15° onwards, 
it increases rapidly to ~ ID. At the same angle 9 — 
15° = 9d, we observe a very spectacular transition in the 
evolution of £™ ax . The latter first increases continuously 
from 51? at 9 — to 15D at 9 = 9d- Here, a sudden jump 
brings up £™ ax from 15D to 30L>. The cutoff at 30D is 
imposed by data processing for the evaluation of £ c for 
individual grains. However, this value is close enough to 
the depth H = A0D of the granular bed to be interpreted 
as corresponding to a system size correlation. 

It is interesting to map the growth of £™ ax and its 
nearly discontinuous change at 9 = 9d into a percolation 
process of the overloaded grains. Indeed, £™ ax increases 
with 9 because of the growing number of overloaded 
grains (Figure EJ) ■ Each overloaded grain represents an 
active site which can belong to a nearest-neighbour clus- 
ter. In this picture, £™ ax corresponds to the size of the 
largest cluster of active sites, while £™ ean represents the 
mean cluster size. In this sense, the transition at 9 = 9 4 
is a percolation transition where the largest cluster size 
diverges (spanning the whole system). The fraction of 
active sites at the transition is p = 0.45 (Figure |5J|. In- 
terestingly, this proportion is the percolation threshold 
for randomly built-up structures of non-overlapping par- 




FIG. 11: Successive circular neighbourhoods of diameter £ 
centered on one overloaded grain (in gray). The stress ratios 
7f of the grain are computed over these neighbourhoods as a 
function of I. 

tides in two dimensions [23, [3~1, |3^| . 

FigureHHl shows successive snapshots of the packing in 
the frame attached to the simulation box for an increas- 
ing slope angle with overloaded grains filled in black. We 
clearly observe the clusters of overloaded sites that grow 
in size and coalesce already at 9 — 15° to form a con- 
nex network spanning the packing from top to bottom. 
Beyond this point, the overloaded grains continue to per- 
colate (while their number increases) until slope failure 
occurs. 

Note that the circular shape of our control volume 
does not differentiate space directions. However, as men- 
tionned before, because of our tight grain size distribu- 
tions, there are privileged directions at which a larger 
concentration of overloaded grains can be observed (Fig- 
ure El). This means that, if the percolation of over- 
loaded grains was analyzed for different space directions, 
we would get a "directed" percolation at those privileged 
directions for a slope angle slightly below 9 = 15°. 



VI. TRANSITION TO COHERENT SHEARING 

The evolution of the packing stress ratio T with 9 (Fig- 
ure |5j) carries no apparent signature of a discontinuous 
transition at 9d- Irrespective of boundary conditions and 
for different values of the aspect ratio H/L, the packing 
stress follows closely the loading direction up to slope fail- 
ure. Such a signature should thus appear in the volume- 
change behaviour of the packing. Indeed, numerous small 
rearrangements in the bed lead to small variations of the 
total volume V of the bed. Figure ^] shows the volu- 
metric strain sy = {V — Vq)/Vq, where Vq is the initial 
volume, as a function of 9 in the PBC system. We see 
that the volume first decreases (negative values of ey) 
with 9. Then, precisely at 9 = 9d, this contractant be- 
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FIG. 12: Evolution of £™ ean (a) and C™ (b) as a function 
of 9 (see text for definitions). 



haviour transforms into a dilatant behaviour, with the 
volume increasing (positive values of £y) up to surface 
failure. 

Since packing dilatation may occur only as a conse- 
quence of shearing, the transition to dilatation at 64 can 
be interpreted as the onset of a stable shear mode in the 
packing. Before transition, the particle rearrangements 
occur in a diffuse and incoherent manner in the whole 
packing. This is due to the fact that, in the presence 
of geometrical disorder and gravity, each grain tends to 
occupy a more stable position in the packing. Such re- 
arrangments take place mostly in a collective way but in 
small volumes compared to the size of the system [l3j . 
and they do not disturb the overall stability of the slope 
(neither at the free surface nor in the bulk) as long as 
they have low amplitudes (involving low rotation rates). 

This behaviour is reminiscent of the contractant be- 
haviour observed in the first stages of a shear test per- 
formed on soil samples |4jj • The extent of volume reduc- 
tion depends on the initial solid fraction of the packing. 
A transition to dilatancy happens if the initial solid frac- 
tion is above the "critical state" solid fraction, i.e. the 
solid fraction corresponding to a state reached after long 
enough monotonous shearing 3]. This is an interesting 
analogy although in a granular slope the gravity behaves 
as a bulk force in contrast to shear testing conditions 
where the largest stresses are exerted at the wall bound- 
aries on the sample. Moreover, the rotation of principal 
stress directions with respect to the packing (as in our 
rotating bed) is not a "standard" testing method for the 
behaviour of soil samples. 




FIG. 13: Maps of the overloaded grains (in black) for suc- 
cessive values of 9. From top to bottom, 9 — 0, 9 — 5°, 
9 = 10°, 9 — 15° and 9 — 9 C . They form growing clusters 
which eventually span the packing from top to bottom. 




6 (deg) 

FIG. 14: Evolution of the volumetric strain ev = (V — Vo)/Vo, 
Vo being the initial volume, as a function of 9 in the PBC 
system. 



VII. DISCUSSION AND CONCLUSION 

A. Summary 

In this paper, we analyzed the stress states in a gran- 
ular slope composed of rigid disks and simulated by the 
Contact Dynamics method. We showed that a global 
stress analysis in terms of Euler's equations for periodic 
boundary conditions along the free surface fits excellently 
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the simulation data. The influence of side walls was also 
evaluated and compared to the data as a function of the 
total length L of the free surface. 

The local stresses were described in terms of internal 
moments which allow for an additive coarse-graining of 
stresses up to the system size. We observed a wide distri- 
bution of grain stress ratios involving a significant frac- 
tion of the so-called "overloaded" grains, for which the 
stress ratio 7 is above the "critical" stress ratio r c , cor- 
responding to the Coulomb yield threshold for the whole 
packing. 

As long as the slope angle 9 is below the angle of re- 
pose 9 C , the stress ratio in a control volume centered 
on an overloaded grain decreases with the size £ of the 
control volume. This property allowed us to define a 
"correlation" length l c as the size of a control volume 
centered on each overloaded grains, and for which the 
stress ratio is just equal to the critical stress ratio r c . 
The largest length £™ ax in the packing is thus a coarse- 
graining length beyond which all stress ratios are below 
the Coulomb yield threshold. 

We found that £™ ax increases as the bed is slowly tilted 
towards its angle of repose; at an angle Qd below the an- 
gle of repose t™ ax undergoes a sudden jump to a length 
comparable to the depth of the bed. We argued that this 
evolution of local stresses can be mapped into a perco- 
lation process of the overloaded grains. In this picture, 
the discontinuous increase of the coarse-graining length 
at corresponds to a percolation transition. Remark- 
ably, this transition coincides with the onset of dilatation 
in the packing. 



B. Discussion 

The above results reveal an unexpectedly rich scaling 
of stresses in a granular bed in the vicinity of slope fail- 
ure. In particular, they are consistent with the advent of 
a "phase transition" at 9d- Upon transition, the stress 
inhomogeneities, in the form of clustering of overloaded 
grains, extend to the whole packing. At the same time, 
the grain rearrangements turn from an incoherent be- 
haviour into coherent dilatant shearing. This transition 
is all the more interesting that it occurs far enough from 
the angle of repose 9 C to be considered as a strong pre- 
cursor of the incoming surface instability. 

Although the transition angle 9d appears here in the 
course of a "quasi-static" evolution of the slope, it is still 
tempting to identify 9 a with the dynamic angle of repose, 
i.e. the angle reached when the slope comes to rest after 
failure. This is an appealing interpretation in that it 
relates the dynamic angle of repose to a static property of 
the packing. In other words, the order parameter in this 
description is the fraction of overloaded grains. Further 
simulations are necessary to check this conjecture, but 
let us simply consider here an argument in this direction. 

In the presence of long-range correlations or large clus- 
ters of overloaded grains for the angles 9 in the interval 



[6d, 9 C ], the granular bed should be metastable (sensitive 
to perturbations). This is consistent with experimental 
measurements of the size of the perturbation (number 
of grains added at a point on top of the bed) required 
to trigger an avalanche, which is shown to increase as 
a function of the distance 59 = 9 C — 9 to the angle of 
repose A surface flow triggered at 9 C acts as a per- 
turbation that will continue to exist as long as the slope is 
metastable. Hence, it will not stop until 9 declines below 
9d- This means that the avalanche destroys long-range 
correlations and takes the slope back to a "disorganized" 
state. 

It is worth noting that we did not observe such a clear 
transition for the pair correlation functions of contact 
forces and grain stresses at 9d- In these cases, the corre- 
lation length remains of the order of a few grain diame- 
ters all along the simulation. However, as shown in detail 
in , the evolution of "critical contacts", i.e. contacts 
where the friction is fully mobilized, is consistent with 
a percolation- like process (of critical contacts), leading 
to a power-law increase of the corresponding correlation 
length from 9d onward and a divergence of the correlation 
length at 9 C . This behaviour is plausible in view of the 
onset of coherent shearing at 9d, as shown in this paper. 

This suggests that the evolution of a granular bed 
is best described in terms of overloaded grains in the 
range [0, 64] and in terms of critical contacts in the range 
[9d,9 c ). Indeed, at slope angles below 9d, the fraction of 
overloaded grains is a good parameter in that it increases 
strongly with the slope angle and it controls the percola- 
tion transition at 9d- In this interval, the friction mobi- 
lization is rather chaotic and incoherent, and the fraction 
of critical contacts increases very slowly with the slope 
angle [13|, At slope angles above 6^, due to coherent 
shearing, the friction forces are increasingly and coher- 
ently mobilized so that the evolution can be described 
in terms of critical contacts (the order parameter being 
the fraction of these contacts). This evolution leads to a 
power-law divergence of the correlation length for critical 
contacts at 9 C . 

The dilatation of the bed provides another interest- 
ing insight into the behaviour of a granular bed in the 
metastable range. The avalanche can not occur unless 
the bed dilates. Since the bed begins to dilate at 9d, the 
angle 9 C — 9d should be interpreted as the "dilatation an- 
gle" of the granular material when plastified as a result of 
rearrangements induced by the rotation of the bed |36| . 
Hence, in a macroscopic approach, the static angle of re- 
pose 9 C and the transition angle 9d, interpreted as the 
dynamic angle of repose, may be described in terms of 
the Coulomb yield criterion and the dilatation angle, re- 
spectively. 
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